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ABSTRACT 

Aims. We present 11 high-precision photometric transit observations of the transiting super-Earth planet GJ 1214b. 
Combining these data with observations from other authors, we investigate the ephemeris for possible signs of transit 
timing variations (TTVs) using a Bayesian approach. 

Methods. The observations were obtained using telescope-defocusing techniques, and achieve a high precision with 
random errors in the photometry as low as 1 mmag per point. To investigate the possibility of TTVs in the light curve, 
we calculate the overall probability of a TTV signal using Bayesian methods. 

Results. The observations are used to determine the photometric parameters and the physical properties of the GJ 1214 
system. Our results are in good agreement with published values. Individual times of mid-transit are measured with 
uncertainties as low as 10 s, allowing us to reduce the uncertainty in the orbital period by a factor of two. 
Conclusions. A Bayesian analysis reveals that it is highly improbable that the observed transit times show a TTV, when 
compared with the simpler alternative of a linear ephemeris. 



Key words, stars: planetary systems, stars: individual: GJ1214, methods: statistical, methods: observational, techniques: 
photometric 



1. Introduction 



2. Observations and Data Reduction 



The transiting exoplanet GJ 1214 b was discovered i n 2009 
by the MEarth projecf] ( |Charbonneau et al.||2008|). This 
plane t transits a nearby M dwarf ( Charbonneau et al. 
20091, with a mass 0.15 Mq and the planet is gener- 
ally classified as a super-Earth with a mass and radius 
(M D = 6.37 M m and R p = 2.74 R e according tolKundurthy 

= 



et al. (|201l|), in this study we find 
77 



M„ 



and 



2.85 R®) between that of Earth and Neptune, a type 
of planet that has no solar system analogue. GJ 1214 b is 
one of the lowest temperature transiting exoplanets known 
and, as it is also detectable with radial velocity (RV) meth- 
ods - a very interesting target for detailed study. 

Due to the relatively low mean density of G J 1214 b, 
p p = 1.49 ± 0.33 gem -3 in this study, it has been suggested 
to hold some extended atmosphere or gaseous envelope. But 
the plane t composition in th is mass and radius range is de- 
generate ( | Adams et al.|200"8 ), warranting further studies in 
order to determine its composition. Defining what the at- 
mosph ere consists of can help d etermine the planet compo- 
sition. Rogers & Seager (20101 describes three possibilities 
for the interior and atmospheric composition of GJ 1214 b. 
It could be (1) a mini-Neptune with a H/He gas envelope, 
(2) a "water world" with a water-rich and ice-dominated 
interior and a water-vapour-dominated envelope, or (3) a 
rocky planet with an atmosphere mainly consisting of H 2 . 

Recent results fr om among others the Kepler mission 
dLatham et al.|20"TT ) and gravitational microlensing ( |Gould 
et al. 2010[ ) gives reason to believe that multiple systems 
are common. It is therefore inherently interesting to look for 
traces of transit timing variations in any transiting system 
and especially so in GJ 1214, as the planet seems to be at 
the inner edge of the habitable zone, with an equilibrium 
temperature of in the region of 393 to 555 K ( Charbonneau 



et al. 2009), i.e. finding a planet in a slightly larger orbit 
would be very interesting. 

Additional planets can be revealed via their gravita- 
tional effects on the transiting planet. This would result in 
telltale systematic deviations in the mid-transit times from 
a linear ephemcris, a phenomenon known as transit tim- 
ing v ariation (TTV) (|Agol et al.||2005| |Holman & Murray 



2005 1. GJ 1214b is well-suited to such an analysis due to its 



short transits, allowing precise measurements of mid-transit 
times, and its short period, which means many transits are 
observable. One disadvantage is its relative faintness, which 
could cause a loss of precision in the measured transit times. 

We have gathered photometric observations of 11 tran- 
sits in the GJ 1214 system, and modelled them to estimate 
the transit times. Inclusion of results from the literature 
leads to a total time interval of 833 days over which TTVs 
can be investigated. In the following work we cast the prob- 
lem of detecting a TTV signal as a model selection problem 
and via Bayesian methods calculate the probability that the 
data, as a whole, actually contain a TTV signal. 

In section 2 we present the new data and their reduc- 
tion, which in section 3 are analysed in order to obtain new 
transit times and physical properties. Section 4 contains a 
description of the Bayesian model selection process. 



* by the MiNDSTEp collaboration from the Danish 1.54m 
telescope at the ESO La Silla Observatory 
** Royal Society University Research Fellow 

1 https : //www. cf a. harvard. edu/~zberta/mearth/ 



We observed five transits of GJ 1214 b in the period of 
2010 July 6th to 2011 June 3rd, using the Danish 1.54 m 
telescope at ESO La Silla and the focal-reducing camera 
DFOSC. The plate scale of DFOSC is 0.39" per pixel. The 
full field of view is 13.7' x 13.7', but for each transit obser- 
vation the CCD was windowed down to reduce the readout 
time from around 90 s to approximately 30 s. The four tran- 
sits in 2010 were observed through a Cousins I filter, while 
the 2011 transit was observed in Johnson R. An observing 
log is given in Table[l] 

The transits were observed with the telescope defo- 
cused, in order to use longer exposure times whilst avoiding 
CCD saturation. This approach allowed us to decrease the 
Poisson and scintillation noise by e xposing for a larger frac - 
tion of the time during transit (see Southworth et al!p0 09 ) . 
The impact of flat-fielding errors was minimised by the use 
of defocusing and by autoguiding the telescope through- 
out the observations. The diameters of the defocused point 
spread functions (PSFs) ranged from 31 to 42 pixels. 

We also observed one transit of GJ 1214 b using the 
1.52 m Cassini Telescope at the Loiano Observatory, Italy. 
The BFOSC CCD imager was used, and was defocused with 
the same ap proach as taken for the D anish telescope obser- 
vations (see Southworth et alT]|2010 ) . 

One transit was observed simultaneously in the g, r 
and z filters using the Calar Alto 2.2 m telescope and the 
BUSCA four-beam CCD imager. A fourth dataset was ob- 
tained in the u-band but, as expected, yielded a light curve 
which was too noisy to be useful. For further discussion 
on the observing strateg y we used for BUSCA please see 
Southworth et al.| ( |2012| . 

Finally, one transit was observed in 2010 using the 
GROND seven-beam imager on the 2.2 m MPI telescope 
at ESO La Silla. The observations were performed without 
telescope defocusing. Useful results could be obtained for 
only the r and i filters, with the star proving too faint in 
the g, z and JHK channels. 

The data were redu ced using the pipeline described in 
mthworth et al. ( |2009 ) , which performs standard aperture 
photometry with the ASTROLIB / APEeJ^Jidl routine. A range 
of aperture sizes were tried, and the ones which gave the 
least noisy light curves were adopted (Table[l]). The form of 
the transit is insensitive to the choice of aperture sizes, and 
to the presence of two faint stars within the sky annulus. 

Comparison stars to be used for relative photometry 
were chosen within the field, and the ASTROLIb/aper rou- 
tine determined relative magnitudes for these and GJ 1214, 
from the given coordinates and aperture radius. Potential 
comparison stars that proved to be variable or too faint 
were discarded. Relative photometry of GJ 1214 was ob- 
tained against an optimally weighted ensemble of compar- 
ison stars. 

The resulting GJ 1214 light curve does not have a con- 
stant magnitude out of transit, primarily due to changes 
in airmass and intrinsic stellar variability. To correct any 
systematic trends, the out-of-transit data points were fit- 
ted with a straight line. Simultaneous optimisation of the 
comparison star weights and the out-of-transit polynomial 
was used to obtain the final light curves, which are shown 
in Fig.[l] 



Welcome .html 



2 Distributed within the idl Astronomy User's Library at 
http : // idlastro . gsf c .nasa. gov/ 
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Table 1: Log of the transit observations of GJ 1214 for this work. 



Date 


Telescope 


Start 


End 


Number of 


Exposure 


Filter 


Airmass 


bcatter 


Aperture 


PSF area 




instrument 


(UTC) 


(UTC) 


exposures 


time (s) 






mmag 


sizes™ (px) 


(PX 2 ) 


2010/07/06 


DFOSC 


05:01 


06:44 


33 


150 


J 


1.40-2.18 


0.97 


21,57,77 


1385 


2010/07/14 


DFOSC 


00:56 


05:59 


142 


80-90 


J 


1.33-2.02 


1.18 


19,30,46 


1134 


2010/08/02 


DFOSC 


00:25 


04:31 


154 


60 


/ 


1.24-1.87 


0.99 


15.5,28,46 


755 


oni n /no /o i 
2010/08/21 


1 11 ( \ O / ' 

JJr UsU 


UU:U2 


no . c o 
03:58 


154 


1-: A 

bU 


1 


1 1 1 O A C 

1. AY— 2. AO 


i on 
1.29 


16,26,48 


on a 
804 


2011/06/03 


DFOSC 


02:33 


05:11 


101 


60 


R 


1.61-1.21 


1.36 


16,24,45 


804 


2010/08/06 


BFOSC 


19:52 


22:55 


71 


90 


i 


1.30-1.61 


1.74 


12,20,40 


452 


2011/08/25 


BUSCA 


19:52 


22:58 


123 


90 


9 


1.19-1.30 


3.20 


8,16,24 


201 


2011/08/25 


BUSCA 


19:52 


22:58 


125 


90 


r 


1.19-1.30 


1.92 


20,55,65 


1257 


2011/08/25 


BUSCA 


19:52 


22:58 


125 


90 


z 


1.19-1.30 


2.66 


25,35,45 


1963 


2010/04/29 


GROND 


05:52 


08:38 


89 


50 


i 


1.31-2.54 


3.65 


no defocus 




2010/04/29 


GROND 


05:52 


08:38 


90 


50 


r 


1.31-2.54 


4.80 


no defocus 





Notes. The three numbers are the apertures radii in pixels of the object aperture and the inner and outer edge of the sky 
annulus. 



Three additional light cu rve were obtained fr om the 
Exoplanet Transit Databas^] ( |Poddany et al.|2010 ). These 
were contributed by Johannes Ohiert (2010/07/18) and 
Thomas Sauer (2010/06/29 and 2010/07/07). 



3. Light curve analysis 

3.1. Analysis with jktebop 

The light curves were anal ysed with the jktebof]^] code 
(Sout hworth et al. 2004a b ), originall y developed as ebop 
Popper & Etze^ 



1981 



curves of detached 



Etzel 19811 for modelling light 



eclipsing binaries. The use of the code 



for exoplanet tran sit light curves is discussed in detail in 
Southworth|(|2008|. 



The size of the planet relative to the size of the star 
is directly related to the transit depth, jktebop models 
the sky projections of the two objects as biaxial spheroids, 
dividing them into concentric circles and assigning a limb 
darkening to each of the rings before estimating the flux. 
With optical observations of a planetary system, it is safe to 
assume that the secondary object, the exoplanet, is dark, so 
the surface brightness ratio can thus be set to zero. Using 
jktebop we fitted for the inclination i of the orbit, the 
sum and ratio of the fractional radii k = r p + r± and rp/r*. 
The fractional radii are r p = R p /a and r* = R+/a, where 

R p and i?* are the absolute planetary and stellar radii, 3.2. Physical properties of the system 



2004) provides LDCs for stars as cool as GJ1214A. We 
used values for a star of T c ff — 3000 K and log g — 5.00 (cgs 
units). 

We analysed the combined 2010 Danish telescope data, 
finding that one LD coefficient could be included as a fit- 
ted parameter. We therefore fitted for the linear coefficient 
whilst holding the nonlinear coefficient fixed at the the- 
oretical value. For the other datasets we had to fix both 
coefficients in order to avoid unphysical results. The best 
fits to the Danish telescope data are plotted in Fig. [2] 

In order to obtain uncertainties in the fitted parameters 
we first rescaled the error bars of the datapoints to give 
a reduced x 2 of unity for each transit light curve. This 
step is necessary because the data errors returned by the 
APER algorithm are normally too s mall. We then performed 
1000 Monte Carlo simulations (see Southworth et al. 2004b ) 
on each light curve to derive the errorbars in the fitted 
parameters quoted in Table[2] 

The light curve from 2010/07/06 was unintentionally 
obtained using an exposure time of 150 s, which is signif- 
icant compared to the duration of the ingress and egress. 
When fitting these data we used the possibility within jk- 
tebop to numerically integrat e the light curve m odel in 
order to obtain an unbiased fit ( Southworth||201 1 ) . 



respectively, and a is the semi-major axis of the orbit. We 
also fitted for the time of minimum of each light curve, T m j^ 
using a fixe d orbital period of P = 1.58040490 days (Berta 



et al. 2011). The mass ratio, which only affects the shape 
of the ellipsoids describing the components, was fixed at 
q = 0.0002, which is sufficiently close to the value 0.00013 
found in this paper. We found that changes, less than one 
order of magnitude, in this value have a negligible effect on 
the results. 

Limb darkening (LD) affects the shape of a transit light 
curve. LD will cause the bottom of the transit to have a 
curved shape. We tried fitting four different LD laws: lin- 
ear, square-root, logarithmic and quadratic. LD coefficients 
can be found from stellar atmosphere model s given t he ef- 
fective temperature and surface gravity. Only Claret ( 2000 



It is possible to calculate the physical properties of the 
GJ 1214 system from our measured photometric parame- 
ters and from published values of the stellar effective tem- 
perature, metal abundance and orbital velocity amplitude. 
We obtained final values for r*, r p and i from our results 
for the 2010 Danish Telescope data. Each was calculated as 
the mean of the values from the solutions using the four LD 
laws, with uncertainty the quadrature sum of the largest 
individual uncertainty plus the standard deviation of the 
four individual parameter values. For the star we adopted 
the temperature T c g = 3026 ± 150 K and orbital velocity 
ampli tude K+ — 12.2 ± 1.6ms" 1 from Charbonneau et al. 
(|2009[), and the metal abu ndance [§] = +0.39 ±0.15 from 
Kojas-Ayaia et al.l (|2010|). 



3 http : //var2 . astro . cz/ETD/index . php 

4 jktebop is written in FORTRAN77 and the source code is 
available at |http : / /www . astro . keele . ac . uk/~ jkt/| 



The physical properties were then calculated by requir- 
ing the properties of the star to match the tabulat ed pre- 
dictions of the DSEP stellar evolutionary models (Dotter 
et al.||2008[ ). This step was performed using the procedure 
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Table 2: Photometric parameters of GJ 1214 from the best-fitting light curves to the 2010-season data from the Danish 
telescope, a is the rms scatter of the data around the best fit. 



T T \ 1 , . , , . 

LL> law 


Linear 


bquare-root 


Logarithmic 


Quadratic 


r* + r v 


0.0812 ±0.0030 


0.0803 ± 0.0032 


0.0805 ± 0.0033 


0.0777 ± 0.0040 


k 
i 

u 


0.1222 ±0.0014 
87.97«; 3 2 ^ 
0.45 ± 0.06 


0.1209 ±0.0015 

»S 1 O+0.42 
oo.1Z_q 22 

0.05 ±0^06 


0.1212 ±0.0016 
B8.09±g-g 
0.59 ±0.06 


0.1190 ± 0.0021 
88.501°;™ 
0.28 ±0.07 


V 




0.70 fixed 


0.20 fixed 


0.40 fixed 


r* 


0.0724 ± 0.0027 


0.0717 ±0.0028 


0.0718 ± 0.0028 


0.0694 ± 0.0035 


r v 


0.00886 ± 0.00039 


0.00866 ± 0.00042 


0.00870 ± 0.00044 


0.00826 ± 0.00055 


a (mmag) 


1.18 


1.18 


1.18 


1.18 



Table 3: Derived physical properties of the G J 1214 system. 
The upper part of the table contains the input parameters 
to the property-calculation algorithm and the lower part 
the output parameters. 



' p 

i (deg) 
Teff (K) 

[f ] (dex) 
K+ (ms" 1 ) 



0.0713 
0.00862 
88.17 
3026 
+0.39 
12.2 



±0.0037 

±0.00059 

±0.54 

±150 

±0.15 

±1.6 



M* (M Q ) 
R* (Rq) 
logs* (cgs) 
M p (M Jup ) 
Rp (Rjup) 
9p (ms" 2 ) 

pp (p.Tup) 

a (AU) 



0.150 
0.216 
4.944 
0.0197 
0.254 
7.6 
1.12 
0.01411 



±0.011 

±0.012 

±0.013 

±0.0027 

±0.018 

±1.5 

±0.25 

±0.00032 



outlined by |Southworth| ( |2009| ). The DSEP model set was 
chosen because it is the only one of the five sets used by 
Southworth (2010) which reaches to sufficiently low stellar 



masses. 



The input and output parameters for this analysis are 
given in Table[3] and show a reasonable agreement with lit- 
erature values. The mass, radius and surface gravity of the 
star are given by -M*, i?* and log g+, respectively. The mass, 
radius, surface gravity and mean density of the planet are 
denoted by M p , i? p , g p and p pi respectively. 

The measured physical properties will be subject to sys- 
tematic errors as theoretical evo lutionary models a re not 
perfect representations of reality. Southworth ( 2009 1 found 
that this systematic error was generally 1% or less for the 
masses and radii of transiting planets and their host stars. 
In the case of GJ 1214 this systematic error could be sig- 
nificantly larger, due to the relatively poorer theoretical 
understanding of 0.2 M^stars, but will still be significantly 
smaller than the statistical uncertainties quoted in Table [3] 

3.3. Orbital period determination 

The times are given in Barycentric Julian days (BJD), and 
have been calculated from UTC with codes provided by 
Eastman et al. (20101. By augmenting our measured T mi d 
values with ones from the literature (Table [2]), we are able 
to refine the orbital ephemeris for GJ 1214, refitting for P 
and Tq the zero epoch. Taking the second observed transit 



to be the reference epoch we find the ephemeris: 
T m id = BJD(TDB) 2455320.535733 ± 2.1 • 10~ 5 
+ 1.58040456 ± 1.6 • 10~ 7 x E 



(1) 



where E is the orbit count with respect to the reference 
epoch. The reduced x 2 f° r this fit is 1.24. The estimated 
period is identical to the previous estimates, but the un- 
certainty on the period P has been reduced by approx 
imately a factor of two (|Berta et al 



2011 



20111 Kundurthy et al.l I2011D. The fit and the residuals 



Carter et al. 



are plotted in Fig. [3] Given the rather large spread in data 
points compared with the period, one would not expect this 
estimate of the period to be afflicted by significant system- 
atic error. On the other hand the T could have systematic 
error, which is handled in the later Bayesian analysis by 
marginalising out this parameter. 

4. Transit timing variation as Bayesian model 
selection 

A system with a transiting planet enables the possibility of 
detecting additional unseen planets in the timing data by 
TTVs, that is, look for systematic trends in the residuals of 
Fig. [3] We will in the following outline a method for quan- 
titatively assessing the probability of whether such a signal 
exists in the timing dataset or not. 

One could simply fit an appropriate function describing 
mutual gravitational perturbations, and evaluate a measure 
like the classical squared sum of residuals. But it is, in gen- 
eral, true that a model with more free parameters will be 
able to fit noise better, i.e. produce a lower squared sum 
of residuals by fitting nonphysical features. This is the con- 
cept of over-fitting. The problem is especially prominent 
when the sought effect is on the same order as the uncer- 
tainty in the data, in which case it is not easy to determine 
how much of an improvement in the squared sum of resid- 
uals one should demand for a more complex model to be 
plausible. 

One is forced to penalise models for having too many 
parameters. Thus the problem is no longer a problem of pa- 
rameter estimation, but a problem of model selection. One 
needs to apply Occam's razor. One formal way of doing so 
is via Bayesian model selection, which is completely analo- 
gous to, but distinct from, Bayesian parameter estimation. 
Whilst it might not be possible to estimate the parameters 
of a model with any certainty, a wide range of plausible pa- 
rameters which do improve the squared residuals by some 
amount would lend credibility to the notion that the model 
in question is true. How true can be expressed as a proba- 
bility. 
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Fig. 1: The 11 observed transit light curves of GJ 1214 b, 
plotted in the same order as listed in the observing log in 
Table [T] For the observations taken in multiple filters si- 
multaneously the datapoints are coloured in spectral order, 
i.e. blue for bluest etc. 



The search for a TTV signal can conveniently be cast as 
a model selection problem. If there is a TTV signal there 
has to be some kind of pattern in the residuals listed in 
Table[4] The expected signal from perturbation from an- 
other body can be calculated with an N-body orbital me- 
chanics code. 
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Fig. 2: Plot of the combined Danish telescope light curve 
of the 2010 data (top) and the 2011 light curve bottom), 
versus the best jktebop fit (solid lines 2010, dashed 2011) 
with the logarithmic LD law. The residuals to the fits are 
plotted below, where the lines mark zero residual. 



4.1. Bayesian estimation 
4.1.1. Parameter estimation 



Following the development in (Chapter 3 Gregory 2005), 
Bayesian parameter estimation can conceptually Be 
thought of as testing a range of mutually excluding hy- 
potheses Hi, i.e. one assumes a parametrised model M. 
This model can be thought of as a logical disjunction ("or") 
M = Hi + H% + ■ ■ ■ where the hypothesis Hi implies that 
a given parameter 8 has the particular value 9i. 

For mutually exclusive probabilities the sum rule ap- 
plies. Assuming that the parameter 6 does indeed take a 
value, one can write 



J2p(Hi\M) = l 



From Bayes' theorem one learns that 

p(Hi\M)p(D\Hi, M) 



p(Hi\D,M) 



p(D\M) 



(2) 



(3) 



where the left hand side p(Hi\D, M) is the posterior prob- 
ability of Hi, i.e. the probability of the hypothesis Hi in 
the light of the data D. The term p(D\Hi, M) is the prob- 
ability of the data D if Hi true. This quantity is known as 
the likelihood of Hi. The term p(Hi\M) is known as the 
prior and represents the probability one assigns to Hi in 
the light of one's model before any data becomes available. 
Finally the term in the denominator p(D\M) is known as 
the global likelihood or the evidence. 

Based on assumptions one can deduce the value of 
p(D\M) 



^p{ Hi \D,M) = 



= 1 



Thus 



p(D\M) = ^p(Hi\M)p(D\Hi,M) 



(4) 



(5) 
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Table 4: Mid-transit times from the literature and the 



presen t w ork. Reference 1: T his wor k, 2: Kundurthy et al. 



(2011), 3: Berta et al. (20111 and 4: Carter et al. (2011) 





Epoch Residual (O-C) Ref. 


BJD(TDB) 


xl(T 4 



2454980 


748682 


± 





000104 


-215 


-0 


702161 


3 


2454980 


748570 


± 





000150 


-215 


-1 


822165 


4 


2454983 


909820 


± 


o 


000160 


-213 


2 


586553 


4 


2454983 


909507 


± 





000090 


-213 


-0 


543450 


3 


2454988 


650808 


± 


o 


000049 


-210 





329618 


4 


2454999 


713448 


± 





000115 


-203 


-1 


589879 


3 


2455002 


874670 


± 





000190 


-201 


2 


538827 


4 


2455269 


962990 


± 





000160 


-32 


2 


025110 


4 


2455288 


928200 


± 





001100 


-20 


5 


577393 


4 


2455296 


830130 


± 





000230 


-15 


4 


649176 


4 


2455307 


892454 


± 





000271 


-8 


-0 


430327 


2 


2455315 


794693 


± 





000080 


-3 


1 


731454 


4 


2455315 


794850 


± 





000230 


-3 


3 


301455 


4 


2455315 


794968 


± 





000930 


-3 


4 


481454 


1 


2455315 


795050 


± 





000660 


-3 


5 


301451 


1 


2455315 


794564 


± 





000066 


-3 





441452 


3 


2455318 


955230 


± 





000170 


-1 


-0 


989833 


4 


2455326 


857404 


± 





000110 


4 





521950 


3 


2455334 


759334 


± 





000066 


9 


-0 


406266 


3 


2455353 


724539 


± 


o 


000307 


21 


3 


096014 


2 


2455353 


723870 


± 





000180 


21 


-3 


593988 


4 


2455356 


884950 


± 


o 


000150 


23 


-0 


885273 


4 


2455364 


787000 


± 





000150 


28 


-0 


613490 


4 


2455375 


849970 


± 


o 


000130 


35 





767009 


4 


2455377 


431461 


± 





000420 


36 


11 


631362 


1 


2455383 


752050 


± 


o 


000130 


40 


1 


338790 


4 


2455383 


751635 


± 





000160 


40 


-2 


811207 


1 


2455383 


752143 


± 





000260 


40 


2 


268790 


2 


2455385 


332107 


± 





000610 


41 


-2 


136850 


1 


2455391 


654105 


± 





000059 


45 


1 


660576 


4 


2455391 


654029 


± 





000160 


45 





900575 


1 


2455396 


395438 


± 





000200 


48 


2 


853647 


1 


2455410 


618895 


± 





000100 


57 


1 


012855 


1 


2455415 


360098 


± 





000180 


60 





905925 


1 


2455429 


583692 


± 





000160 


69 





435133 


1 


2455715 


637033 


± 





000170 


250 


1 


583700 


1 


2455799 


397830 


± 





000500 


303 


-4 


865397 


1 


2455799 


398465 


± 





000270 


303 


1 


484603 


1 


2455799 


398686 


± 





000270 


303 


3 


694603 


1 



That is, the denominator, which does not depend on the 
individual hypotheses Hi, is the sum of the numerator 
over Hi. The process of summing over all the hypothesis 
is known as marginalisation. It is clear that the evidence 
serves as a normalisation constant. For parameter estima- 
tion this normalisation is unimportant, as in most cases one 
simply seeks the maximum posterior value without regard 
to normalisation. But for model selection this evidence term 
is important. 



4.1.2. Model selection 

Model selection is carried out following the exact same 
procedure as above, only assuming a disjunction of mod- 
els / = Mi + M2 + ■ ■ ■ + MjVj instead of a disjunction of 
hypotheses. 

Bayes' theorem now reads 



p{Mi\D,I) 



p(M\I)p(D\M„T) 
p(D\I) 



(6) 



One immediately recognises the second term in the numera- 
tor as the evidence term from Eq.[5j Just as the probability 
of a parameter is proportional to the likelihood times the 
prior, the probability of a model is proportional to the ev- 
idence times the prior. The denominator in Eq.[6] can be 
calculated like the evidence in Eq.[5j assuming that one of 
the Mi is the true model. In other words the probability of 
a model is given by marginalising over all the parameters 
in the model. 

Often it is more convenient to work with odds ratios O 
and Bayes factors B defined as 



O, 



p(Mj\I) p(D\Mj,I) 
piM^piDlMjJ) 



P(Mj\I) 
V(M 3 \I) ij 



Assuming that 



(7) 



(8) 



N 



one can write the probabilities of the individual models as 



p(Mi\D,I) = 



J2n °ji 



(9) 



It is crucial to notice that these odds ratios implicitly re- 
sult in an Occam's razor. One may notice that the evidence 
term in Eq.[5] takes the form of an average of the likeli- 
hood over the prior. Hence penalising complicated models 
by their unused prior space, e.g. large swathes of prior space 
with negligible likelihood, will pull down the average of the 
likelihood over the prior. In other words, the more plau- 
sible model is the one that makes more sharp predictions 
with fewer adjustable parameters. Note that the form of 
the priors and the prior ranges are as much part of the 
model specification as the functional form of the likelihood 
function, which is no surprise given that probabilities are 
purely a function of one's state of knowledge. There is noth- 
ing inherent in nature about probabilities. Also note that 
assigning probabilities according to Eq.[9j i.e. normalising 
to a total sum of one, implicitly assumes that one of the 
models is true. 



4.1.3. MultiNest 

The mainstay of Bayesian posterior distribution evaluation 
has for many years been Monte Carlo methods, in par- 
ticular Monte Carlo Markov Chain (MCMC) algorithms. 
Unfortunately it turns out that most MCMC incarnations 
generally have trouble accurately estimating the evidence 
term, in particular when the posterior has multiple modes. 
Indeed many of the various MCMC packages available do 
not directly calculate the evidence term. 

As noted before, the estimation of parameters does not 
require one to calculate the evidence term, but model selec- 
tion does. A recent innovation in Monte Carlo methods is 
the MultiNest algorithm, whi ch is specifically designed to 
accurately evaluate evidences ( Feroz fc Hobson|2008 Feroz 
et~aI7|[2009l ) . 



4.2. Application of Bayesian model selection to TTV search 

TTVs arise from the dynamics of a planetary system if 
there is more than one planet in the system. But as it is 
commonly known, the three body problem and higher is 
chaotic over long time scales. Such a signal could show up 
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Fig. 3: Times of mid-transit from the literature and this work. In the upper plot the errorbars are smaller than the point 
marker and are thus not shown. The lower plot shows the residuals on a different scale. Some of the observations are 
taken with instruments that simultaneously observe in several filters so some of the data points coincide in time. 



as a change in apparent orbital period, or the phase of the 
orbit. 

Given these ambiguities we have then chosen to frame 
the search for a TTV signal as a Bayesian model selection 
problem. To every transit we observe we can unambiguously 
assign an epoch, given prior information on the orbital pe- 
riod, which in the case of GJ 1214 b is approximately 1.5804 
days. We have calculated the probability of the linear no- 
TTV model versus models with TTV simulated with an 
N-body code. 



4.2.1. Non-TTV model 

Assuming no transit timing variation, we would expect that 
the relation between the time of mid-transit and the epoch 
is perfectly linear. In the Bayesian framework presented 
above, given the distribution of errors and the prior ranges, 
we can calculate the probability of this particular model 
itself and compare it directly to the probability of the model 
with a TTV signal. 

As the prior ranges are included in the overall posterior 
probabilities of the models it is necessary to assign ranges 
to these priors, i.e. the posterior probability takes the form 
of an average of the likelihood over the prior; see Eq.[5j 



The parameter Tq in both the models is related to the 
definition of the epoch. It is simply the Julian date of 
the mid-transit of epoch 0. When the epoch is defined, 
this quantity is in principle known, but we cannot deter- 
mine it exactly; hence, we introduce a systematic error. 
In a Bayesian framework we can take this into account by 
marginalising out To- The prior range of To has been chosen 
to correspond to the largest error in the two measurements 
that pertain to epoch in Table|4] 

4.2.2. TTV models 



Unfortunately there is no known analytic solution to the 
general three body problem. Hence, to calculate the TTV 
arising to a third perturbing body in the system one has 
to rely on num erical orbital integration codes. One such 
code is Swift (Levison & Duncan 1994), which has been 
employed here. For our purposes we have adap t ed th e 
FORTRAN code used in |Nesvorny fc Morbidelli| ( |2008[ ), 
with this adapted code we where able to satisfactorily re- 
produce results in that article. 

For the purpose of this investigation it was assumed that 
GJ 1214 b is in a circular orbit, wh ich is likely given the 
results in Charbonneau et al. ( 2009 ) where the eccentricity 
is estimated to be less than 0.27. Further we assume that 
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Relative perturber period 

Fig. 4: Plot of the standard deviation, in units of seconds, of 
the TTV signal of GJ 1214 b calculated for different masses 
of the perturber. The various mean motion resonances are 
clearly visible. The blue shaded area marks perturber orbits 
that would make the orbit of GJ 1214 b unstable. 



GJ 1214 b is perturbed by an planet in a coplanar circular 
orbit. 

Given a trial mass m, period P and mean longitude A of 
the hypothetical perturbing planet, and assuming the time 
of first transit to be t = 0, the provided code calculates the 
time of all future transits taking into account the gravita- 
tional interaction between the planets. In the simple case 
of two coplanar circular orbits the argument of periapsis 
simply determines the angular separation of the two plan- 
ets at the start of the integration. The standard deviation 
of the TTV signal as a function of these parameters can be 
calculated numerically as done in Fig. [3] 

The majority of perturber parameter space will give rise 
to negligible TTV, but from Fig. [4] two regions of perturber 
parameter space which can give rise to significant TTV, 
comparable to the uncertainty of our data, can be identi- 
fied. One region consists of periods in the range from the 2:1 
resonance, at a relative period of 0.5, to the inner edge of the 
instability strip. The other region is from the outer edge of 
the instability strip to the 2:1 resonance, at a relative period 
of 2. In both of these regions, perturbers down to a mass 
of about 0.1 Earth mass will give rise to TTV of a second 
or more. Hence we will use these ranges as prior ranges in 
the model selection. The two regions will be treated as two 
different models to be tested against each other. The upper 
mass range for these two mod els will be set by th e accu - 
racy of radial velocity data from Charbonneau et al. ( 2009 1 , 
where radial velocity data with a accuracy of about lOm/s 
were presented. We estimate that in the outer model, per- 
turbers heavier than 6 Earth masses are excluded and for 
the inner model, perturbers heavier than 5 Earth masses 
are excluded by RV data. The two TTV models are pre- 
sented in Table [5] together with the non-TTV (i.e. linear) 
model, and their priors. 



4.3. Discussion 

The final results of these three models have been sum- 
marised in Table[5j along with the prior ranges of the pa- 
rameters. We conclude that there is a very strong preference 
for the no-TTV model, given the calculated overall proba- 
bilities of the three models, because the relatively low qual- 
ity, in the form of large uncertainties, of the available data 
does not allow us to reach complicated conclusions. 

It should be noted that posterior probability of the 
model includes these prior ranges, hence the probability 
of the model can be thought of as the probability that a 
planet in this part of parameter space explains the tran- 
sit times of GJ 1214 b. But note that the probability of a 
model is normalised by the prior volume, so a model with 
larger prior ranges or more parameters will be less proba- 
ble, unless heavily supported by data. Large prior ranges 
are equivalent to lack of prior knowledge, hence one needs 
strong evidence to accept a complicated model in absence 
of prior knowledge. 

MultiNest is in essence a Monte Carlo code and it does 
produce posterior samples, which can be analysed similarly 
to the samples produced by conventional Monte Carlo codes 
to produce estimates of the posterior probability of the pa- 
rameters and their correlations. This has been done in Figs. 
[5] and [6] for the two TTV cases, respectively. 

The posterior density of the period P in the no- 
TTV show s go od agreement with the estimate produced 
in section 13.31 To evaluate the evidence term we have 
marginalised over a range of offsets to take account for the 
uncertainly in the timing of the zero epoch. 

Similarly the plot of the posterior probabilities in Fig. [5] 
and [6] shows the marginal probabilities of the parameters of 
a perturbing planet, assuming that the data is explained by 
a perturbing planet within the prior ranges. These marginal 
probabilities show preference for a light planet (small m). 
The gaps in the marginal posterior of the period P occurs 
at strong resonances that would give rise to TTV much 
larger than what is observed. 



5. Conclusions 

In this paper we have presented and analysed 11 new tran- 
sit light curves of the exoplanet GJ 1214 b. In addition 
new analyses of three previously published transits are pre- 
sented. With these new data it has been possible to improve 
the previously estimated ephemeris and physical properties 
of G J 1214 and GJ 1214 b. The calculated physical prop- 
erties are consistent with previously published values and 
the uncertainty in the orbital period has been reduced by a 
factor of approximately two - P = 1.58040456 ± 1.6 • 10~ 7 . 

Furthermore a stringent analysis, incorporating avail- 
able prior knowledge in the form of orbital mechanics and 
previous RV investigations, of whether transit timing varia- 
tion can explain the data has been undertaken via Bayesian 
methods. Specifically we have calculated the probability of 
three models for explaining the data, two with a TTV and 
one without, assuming a priori that one of the models is 
true. This analysis has revealed that the models with TTVs 
are highly improbable compared to the simple model as- 
suming no TTV. That is, the given data does not allow us 
to conclude that there is a planet in the mass range 0.1-5 
Earth-masses and the period range 0.76-1.23 or 1.91-3.18 
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Table 5: Table of parameters for the non-TTV and two TTV models that have been analysed. 



Parameter 


Prior range 


PrinT - t^mp 
j. iiui ly uc 


Mean 


St. Dev. 


MAP 


In 1 pvinpnpp I 

1111 C VlU-ClltjC I 


j_ luijauiiiiy 


Linear T mid 


= P ■ E + T 










134.98 


0.9999 


T a 


±0.0002 


Uniform 


2455320.535732 




2455320.535533 






P in days 


±0.0000002 


Uniform 


1.580405 


0.000015 


1.580340 






Interior perturer: T mi d = TTV(a,m 


,X,E)+ To 








122.09 


3- 10" b 


T a 


±0.0002 


Uniform 


-0.000003 


0.0001 


0.0001 






X 


- 2tt 


Uniform 


3.1 


1.2 


1.9 






P in days 


0.76-1.23 


Jeffrey 


0.9 


0.1 


0.9 






m in M & b 


0.0000003-0.0000135 


Uniform 


0.0000009 


0.000001 


0.0000006 






Exterior perturber: T m id = TTV(a, 


m, X, E) ± T 








122.71 


5 • 10" B 


T a 


±0.0002 


Uniform 


0.000002 


0.0001 


0.00004 






X 


- 2vr 


Uniform 


3.1 


1.4 


6.0 






P in days 


1.91-3.18 


Jeffrey 


2.7 


0.3 


2.7 






m in Mq' 


0.0000003-0.000018 


Uniform 


0.0000002 


0.000002 


0.0000005 







Notes. Nu isance parameter marginalised out.' 6 ' These are scale parameters and should have a Jeffrey prior according to 
Gregoryl |2005l). 



Marginal 



A in radians 



Conditional 




Conditional 



S i 

"O 
CO 




'0.6-0.4-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

m in solar massesxur 
Marginal 




l.Op 
0.8 

) 

' 0.6 
} 

3 0.4 
i 0.2 
j 0.0 
j -0.2 
"-0.4 

-°4s 



2.0 2.5 3.0 3.5 4.0 

P in days 



Conditional 



m in solar masses 



2.0 2.5 3.0 3.E 

P in days 



,0.00053 
0.00035 
0.00026 
0.0001 1 
0.00005 

■0.00001 

iO.00123 
0.00082 
0.00061 
0.00025 
0.00012 

■0.00001 



Marginal 




2.0 2.2 2.4 2.0 2.8 
P in days 



Fig. 5: Marginal probability and conditional probabilities for the TTV model assuming an exterior perturber from 
MultiNcst. The parameter T is a nuisance parameter, which is integrated out, hence the marginal and conditional 
probabilities have not been plotted. 



days. To be able to reach such conclusions we would need 
many more consistent data points and/or higher accuracy. 

A planet with a greenhouse warming and albedo sim- 
ilar to Earth at a period of approximately 4.5 days in 
the GJ 1214 system, corresponding to a period relative to 



GJ 1214 b of 3, would have a surface temperature of about 
80° C. Conversely a planet in this orbit with an albedo sim- 
ilar to Venus and a greenhouse warming of would have 
a surface temperature of about 0° C. Hence such a planet 
could be at the inner or outer edge of the habitable zone of 
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Fig. 6: Marginal probability and conditional probabilities for the TTV model assuming an interior perturber from 
MultiNcst. The parameter T is a nuisance parameter, which is integrated out, hence the marginal and conditional 
probabilities have not been plotted. 



G J 1214 b depending on the parameters - for high albedo 
and/or low greenhouse warming the habitable zone overlaps 
with the region of parameter space that can reasonably be 
investigated with a transit timing accuracy on the order on 
10 s; see Fig. |4j Because of the orbital resonance, habitable 
planets down to Mars-mass could potentially be revealed in 
the GJ 1214 b system with the already achieved timing ac- 
curacy. To investigate the habitable zone for planets more 
similar to Earth, a much higher accuracy on the order of 
0.01 s is called for. 
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